Simulation studies of dark energy clustering induced by the formation of dark matter 

halos 
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In this paper, we present a simulation method within the two-component spherical collapse model 
to investigate dark energy perturbations associated with the formation of dark matter halos. The 
realistic mass accretion history of a dark matter halo taking into account its fast and slow growth 
is considered by imposing suitable initial conditions and isotropized virializations for the spherical 
collapse process. The dark energy component is treated as a perfect fluid described by two important 
parameters, the equation of state parameter w and the sound speed c s . Quintessence models with 
w > — 1 are analyzed. We adopt the Newtonian gauge to describe the spacetime which is perturbed 
mainly by the formation of a dark matter halo. It is found that the dark energy density perturbation 
Sde depends on w and c s , and its behavior follows closely the gravitational potential 4? of the dark 
matter halo with Sde ~ — (1 + w)&/c%. For w > — 1, the dark energy perturbation presents a 
clustering behavior with Sde > during the entire formation of the dark matter halo, from linear 
to nonlinear and virialized stages. The value of Sde increases with the increase of the halo mass. 
For a cluster of mass M ~ 1O 15 M , Sde ~ 10 ° within the virialized region for c\ £ [0.5, 1], and it 
can reach Sde — 0(1) with c 2 s — 0.00001. For a scalar-field dark energy model, we find that with 
suitably modeled w and c s , its perturbation behavior associated with the nonlinear formation of 
dark matter halos can well be analyzed using the fluid approach, demonstrating the validity of the 
fluid description for dark energy even considering its perturbation in the stage of nonlinear dark 
matter structure formation. 

PACS numbers: 95.36.+x,98.80.-k 



I. INTRODUCTION 

Since the discovery of the accelerating expansion of 
the universe, tremendous progress has been achieved in 
dark energy studies [l[ [2] [3] [4] [5] f|. Future generations 
of cosmological observations are expected to be able to 
provide us tight constraints on properties of dark energy 
H @ G3 ■ The full realization of their constraining 
power, however, depends on our thorough understand- 
ings on how different dark energy models affect the ex- 
pansion of the universe and the formation and evolution 
of cosmic structures differently [HJ [l3[ . 

For dynamical dark energy models, apart from the 
background field that drives the accelerating expansion 
of the universe, dark energy perturbations exist intrin- 
sically. Their behaviors on large scales and the effects 
on the linear power spectrum of the matter density per- 
turbations and the anisotropy of the cosmic microwave 
background radiation (CMB) have been extensively in- 
vestigated and taken into account in derivin g co smolog- 
ical constraints from CMB observations 14 1 To|[T6jp7j. 
On small scales of the nonlinear structure formation, they 
should also be considered. Therefore a fully consistent 
simulation to study the impacts of dark energy on the 
structure formation should simulate the evolution of the 
dark energy component explicitly together with the mat- 
ter components (dark matter and baryonic matter). The 
additional dark energy component can however compli- 
cate the numerical simulation significantly. Therefore in 
many studies, the dark energy is not included as an in- 
dependent component in simulations. Instead, approxi- 
mations or simplifications suitable to specific models are 



often adopted so that the effects of dark energy pertur- 
bations on the structure formation can be partially taken 
into account [l^[l^[2^[2l|[22]|[23|. In a very recent se- 
ries of papers regarding coupled dark matter-dark energy 
models, f(R) gravity models and scalar-tensor theories, 
the corresponding scalar field com pon ent has been ex- 
plicitly included in simulations QlllldijU^il]. Such 
simulations not only can study the structure formation 
self-consistently, but also can provide the detailed dy- 
namical evolution of the relevant scalar field component, 
which may contain additional information in differenti- 
ating different models. 

To study the effects of dark energy and particularly 
the behavior of dark energy perturbations along with the 
formation of nonlinear structures, the analyses under the 
spherical collapse model are extensively performed. Al- 
though it is a simplified model comparing with the full 
cosmological simulations, it allows us to isolate differ- 
ent effects and to explore a broad parameter space ef- 
ficiently. Studies along this line have shown that the 
simple extrapolation from the amplitude of the dark en- 
ergy perturbation obtained in the linear stage of struc- 
ture formation can grossly overestimate the level of the 
dark energy perturbation in the nonlinear epoch of the 
structure formation [29[ 3(J 31 1. For minimally-coupled 
quintessence dark energy models with <~ 1, the dark 
energy perturbation induced by the formation of cluster- 
scale dark matter halos remains at the level of < 10 ~ 5 , 
in contrast with the amplitude of ~ 10 ~ 2 obtained by 
extrapolating the linear analyses to late stages of dark 
matter halo formation where nonlinear processes actu- 
ally happen^ [H S3 [32] 33] [34] 35] [36] 37] 38] [3j| ■ 
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Such a large difference indicates that the virialization 
of dark matter halos plays an important role in deter- 
mining the behavior of dark energy perturbations in the 
halo region. Previous analyses take into account the viri- 
alization process by adding in a simple prescription into 
the analytical formulation to prevent the halo material 
from collapsing all the way to the center (e.g., [30j [3lj ) . 
Although such a modeling can catch the basics of the 
virialization process, it cannot describe it very realisti- 
cally. Particularly, it cannot handle properly the shell 
crossing that is important for the virialization process. 

To overcome the shortcomings of the artificial treat- 
ment of the virialization process, in this paper, we carry 
out numerical simulations for a set of two-component 
spherical systems. Following Lu et al. [4l[, by specifying 
proper initial conditions, the realistic mass assembling 
history for a dark matter halo is considered in the spheri- 
cal 1-D simulations. The tangential velocity is also taken 
into account, which can affect the final density profile 
of a dark matter halo significantly. With these imple- 
mentations, the formation of dark matter halos can be 
suitably simulated, and the corresponding dark energy 
perturbations can be investigated. In our simulations, 
the dark energy component is regarded as an ideal fluid 
specified by two important parameters, the equation of 
state parameter w and the sound speed parameter c s . 
We investigate the dependence of dark energy perturba- 
tions on various quantities, such as w, c s and the mass of 
the dark matter halo. The validity of the fluid approach 
taking into account the dark energy perturbations is also 
analyzed. 

The rest of the paper is organized as follows. In 
Sec. HH we present the formulations related to the two- 
component spherical system. Sec. IIIII describes the 
methodology of simulations. The results and discussions 
are shown in Sec. IIVI and Sec. E] respectively. 



II. FORMULATIONS 

In this study, we simulate a set of two-component 
spherical systems, where the dark matter and dark en- 
ergy components are noncoupling except gravitational 
interactions. Specifically, we follow the evolution of a 
spherically overdense region of dark matter, and analyze 
the induced dark energy perturbations. Because of the 
small amplitudes, we assume that dark energy pertur- 
bations have no effects on the formation of dark matter 
halos [e.g.,[3l|. The detailed methodology for simulating 
the dark matter halo formation will be described in next 
section. Here we present the formulations used in cal- 
culating dark energy perturbations associated with the 
formation of dark matter halos. 

To consider simulations with shell crossing, it is con- 
venient to adopt the Newtonian gauge to describe the 
perturbed spacetime metric, which has been shown to be 
valid even in the stage of nonlinear halo formation [42| . 



The metric can be written as follows, 

ds 2 = -(1 + 2$)d< 2 + a 2 (l - 2V)dx i dx i , (1) 

where a is the cosmic scale factor, and $ and "J are the 
Newtonian potential and the spatial curvature, respec- 
tively. The flat universe is considered here. In the case 
of spherical symmetry, dxidx 1 = dr 2 + r 2 dil. In Fourier 
space, the Einstein's equations in terms of the comoving 
coordinates then read 

k 2 ^ + 3a 2 H[$ + H<I>] = -A-KGa 2 ^8 Pl , (2) 



k 2 [V - $] = 12^Ga 2 ^(pi +Pi)*i, (3) 

i 

where H = a/a is the Hubble parameter, and the dot 
represents the derivative with respect to time t. The 
sum is over different components with i = DM for dark 
matter and i = DE for dark energy. The quantities pi 
and pi are the background energy density and pressure 
for component i, respectively, Spi = (p; — pi), and is 
the anisotropic stress perturbation for component i [43|. 
For the dark energy component, it in general can have 
the anisotropic stress perturbation which can be signifi- 
cant if it is a relativistic fluid (e.g., (44] [45] 46)). On the 
other hand, cosmological observations can only put very 
weak constraints on it, and thus it may be phenomeno- 
logically sufficient for broad classes of dark energy mod- 
els to regard the dark energy component as a perfect 
fluid [H|. Furthermore, it has been shown in [43] that 
for self-interacting scalar-field dark energy models, the 
anisotropic stress perturbation vanishes and thus they 
can be described as perfect fluids. The equivalence be- 
tween a scalar field model and the perfect-fluid descrip- 
tion is also demonstrated by our simulation studies shown 
here in §IV.D. Therefore in our analyses, we treat both 
the dark matter and the dark energy components as per- 
fect fluids with Cj = 0. We then have $ = In the 
following, \I> will be substituted by $. 

On scales of dark matter halos, we have fc 2< & 3> H(<$> + 
H$). Furthermore, 8pDE <C SpoM- Then Eq. (2) re- 
duces to 

fc 2 $ = -A-KGa 2 5p DM . (4) 

Therefore the potential is determined by the dark matter 
distribution only. This simplifies the calculations signifi- 
cantly. 

For the dark energy component, we define 5de = 
Spde/Pde and 9de = ik ■ ^de with vde being the ve- 
locity of the dark energy fluid. Because of the smallncss 
of their amplitudes, we only analyze the linearized dy- 
namical equations for dark energy perturbations, which 
are given by 




= 0, (5) 
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H(l-3w) 



1 



Ode — 



k 2 



1 



■Sde + * = 



(6) 

where w = pde/pde is the equation-of-state parame- 
ter of the dark energy fluid and c s corresponds to the 
sound speed parameter in the rest frame of the dark en- 
ergy fluid and is defined as c 2 = Spde/^Pde- If the 
perturbation is pure adiabatic, we have c 2 s < for the 
dark energy fluid with w < 0, then the perturbation is 
unstable (e.g., 47|H3)- Therefore from physical consid- 
erations, the perturbation of the dark energy fluid can- 
not be pure adiabatic. Here we regard c 2 as a parameter 
[5i|[50] 51] [52] [53|, and analyze how the behavior of dark 
energy perturbations depends on it. 

The above equations form the bases for simulating 
dark energy perturbations. Specifically, the formation 
of a spherical dark matter halo is followed by numerical 
simulations taking into account the effect of background 
dark energy but without including dark energy perturba- 
tions. We then calculate in Fourier space the potential by 
Eq. (U) , and further dark energy perturbations by Eq. ([S]) 
and Eq.®. Finally, dark energy perturbations in real 
space are computed by inverse Fourier transformations. 

In the context of scalar-field dark energy models, their 
fluid description including the presence of perturbations 
relies on finding suitable correspondences between w and 
c 2 and the scalar field [47| . For a scalar field <fi with 
potential V((f>), we have w = Pde / Pde where 

PDE = l<j> 2 + V(<t>), Pde = \P -V{4>). (7) 

The dynamical evolution of the field is given by 

4> + 3H<j) + V = 0, (8) 

where V = dV/d(j). 

The linear perturbation equation is then 



3HS(f> 



T 6<t> + 5V = 4$^>. 



(9) 



The corresponding density and pressure perturbations in 
the fluid description are given by [47j 



Spde 



■ sv, 



(10) 



and 



5p 



DE 



-<j) 2 <$>+<j)5(t)-5V = 5p DE -25V = 5p DE -2V'l 



The divergence of the velocity Ode corresponds to [4 

k 2 , , 



'DE 



a<j) 



(12) 



We then obtain the sound speed-related parameter c 2 by 



a0 



4 = Spde/Spde = 1 + [3.ff (1 - w 2 ) + w] . (13) 



q It is noted that for a given scalar field model, c 2 and 
'w are related to each other. On the other hand, in the 
general fluid description without concerning a particular 
underlying scalar field model, c 2 and w can be regarded 
as two independent parameters that are used to describe 
various models. 

To test the validity of the fluid description includ- 
ing dark energy perturbations, for a specific scalar field 
model with a double-exponential potential given by 



V(0) = V (e a> 



(14) 



where a = 6 and /3 = 0.1 (e.g., [3i(), we simulate the per- 
turbations both for the scalar field explicitly by Eq. ([9]) 
and with the fluid approach with the corresponding w 
given by Eq. ([7]) and c 2 given by Eq. ([TBI . The results 
are then compared. 



III. SIMULATION ALGORITHM 

As discussed in Sec. II, we ignore the feedback effects 
of dark energy perturbations on the dynamical evolution 
of the dark matter component. Therefore in our studies, 
the formation and evolution of a dark matter halo are 
simulated independently of dark energy perturbations. 
Specifically, the dark matter mass of the simulated halo 
is divided into equal-mass shells. These shells are very 
much analogous to mass particles in N-body simulations, 
and we refer them as shell-particles. Their dynamical mo- 
tions are followed using a one-dimensional code adapted 
from Lu et al. [4l| taking into account the effect of the 
background dark energy component. To compute the 
dark energy perturbations induced by the gravitational 
potential of the dark matter halo, at each time step, from 
the positions of the shell-particles, we construct the spa- 
tial dark matter mass density field on regular radial-mesh 
bins with equal width in the radial coordinate r. This 
spherical dark matter density field is then transformed 
into the Fourier space, and the corresponding gravita- 
tional potential is obtained by Eq For each Fourier 
mode, we calculate the dark energy perturbation with 
Eq. ([5]) and Eq. ©. With the inverse Fourier transfor- 
mation, we then obtain the dark energy perturbation in 
the configuration space. 



A. Simulations for the formation of spherical dark 
matter halos 

Although highly simplified, the spherical collapse 
model (SCM) catches the important aspects of halo for- 
mation, such as giving rise to an important collapse 
threshold that can be used to quantify statistically the 
formation for dark matter halos. On the other hand, 
the simple top-hat SCM cannot describe realistically the 
mass assembly of dark matter halos, which involves the 
early-stage fast accretion and merging and late-stage slow 
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growth (e.g., [54j). Furthermore, the analytical treat- 
ment of SCM cannot model the process of virialization 
of dark matter halos naturally. Lu ct al. |4l[ develop a 
one-dimensional numerical simulation method [55|[56| to 
simulate the formation of a spherical dark matter halo. 
By suitably choosing the initial dark matter mass distri- 
bution, the fast- and slow-accretions of the halo forma- 
tion can be properly taken into account. Here we follow 
this strategy to simulate the formation of dark matter 
halos with the important modification to incorporate the 
effect of the background dark energy. 

In the 1-D simulations, the mass of the dark matter 
component is assigned to equal-mass shells, which are 
analogous to particles in N-body simulations. The accel- 
eration a r of a shell-particle at position r depends on the 
matter content within r. Specifically, we have 



H 2 . GM{r)r 

a r = — —il DE (t)(l + 3w)r - — — 
2 [r l 



1)3/2 



I 2 

(15) 



where the first term represents the contribution from the 
background dark energy with the equation of state pa- 
rameter w and the energy density parameter £Ide{P) in 
unit of the critical density of the universe at time t. Note 
that in general w can be time dependent. The second 
term is the gravity from the dark matter component, and 
the third term is the centrifugal force from the angular 
momentum J. Both the Hubble parameter H and the 
dark energy density parameter VLde are time evolving 
and are dependent on the equation of state parameter 
w. In the second term, M(r) = z^2r'<r mr ' ^ s ^ ne total 
mass within r calculated by summing over the mass of 
the shell-particles inside r. The softening parameter a s 
is taken to be a s = 1 kpc, which is about 0.0005i?„i r 
for the virial radius R V i r ~ 2Mpc for clusters of galaxies 
[4lj . The angular momentum J is added in by hand to 
prevent an over-concentrated dark matter density pro- 
file. Following Lu et al. [|l[, a tangential velocity for a 
shell-particle is added in when it falls back to one half of 
its turn-around radius Rt with the tangential and radial 
velocity dispersions, a 2 and a 2 , satisfying the relation 



(16) 



l + (Rt/r a y 



where /3 is taken to be (3 = 2 following Lu et al. [44j 
who show that the results are insensitive to the specific 
value of (3. The parameter r a is chosen to be the virial 
radius of the halo at time a c , the transition time between 
the fast-accretion and slow-accretion phases [see Eq. ([T7| ] 
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Specifically, the radial and tangential velocities are 
randomly generated according to Gaussian distributions 
with the dispersions of a T and at , respectively. Then the 
total kinetic energy of the shell-particle is partitioned 
into radial and tangential components according to the 
square of the ratio of the two random numbers. Note that 
we include the tangential velocity merely in Eq. (|15[) to 
modify the radial motion of a shell-particle without really 
simulating its tangential motion. 



In our simulations, we use 10 5 dark matter shell- 
particles. Further increasing the number does not change 
the results significantly. The symplectic integrator is em- 
ployed ■ The time step is chosen to be smaller than 
the minimum of the dynamical time scales of the shell- 
particles. 

For dark energy perturbations, we do not calculate 
them at the positions of the moving dark matter shell- 
particles. Instead, they are done on a fixed regular radial- 
mesh of equal width along the radial coordinate r in ac- 
cord with the Newtonian-gauge metric shown in Eq. ([1} . 
At each time step of the simulation, we construct the 
spherical dark matter density field on that radial-mesh 
from the positions of dark matter shell-particles. Then 
the over-density field SpoM is obtained and transformed 
into the Fourier space. The gravitational potential $(fe) 
is calculated and the dark energy perturbation for each k 
mode is further solved with Eq. © and Eq. ©. Finally, 
the dark energy perturbation at a desired output time in 
the configuration space is obtained by the inverse Fourier 
transformation. 



B. Initial conditions and virialization 

Apart from the necessity to introduce tangential mo- 
tions, cosmological simulation studies have shown that 
the universal density profile for dark matter halos also de- 
pends sensitively on their mass assembling history, which 
cannot be described well by the simple top-hat spherical 
collapse model [e.g., [5^ . However, as demonstrated by 
Lu et al. [4lj . this can be remedied within the spheri- 
cal collapse framework by setting up properly the initial 
mass distribution for the dark matter component (instead 
of the simple top-hat density profile) according to an ap- 
proximate description about the realistic mass accretion 
process. 

Considering only the growth of the amount of mass for 
a dark matter halo without concerning the detailed pro- 
cesses, Wechsler et al [58| present an approximate form 
for the time-dependence of its mass M(a), which is given 
by 



M(o) = M exp 



-a c S 



(7(1 



1 



(17) 



where Mq is the halo mass at the observed epoch ao, and 
a c is the characteristic scale factor to divide the fast- and 
slow-accretions specified by (din M/dlna)\ a =a c = aoS 
and S = 2. We can see that a c plays the critical role 
in describing the mass growth for a halo, which in turn 
affects the final density profile, such as the concentration 
parameter of the NFW halos [59[[6(|. In other words, 
given a desired concentration parameter for the final den- 
sity profile of a dark matter halo, one can find a suitable 
value of a c so that Eq.(fT7|) can describe the mass assem- 
bly of the halo properly [41| . 

In ACDM models, the empirical relation between the 
concentration parameter and the mass of a halo has been 
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investigated extensively from numerical simulations (e.g., 
[6]j[58{[53|(63|(63i). For example, [H} presents a simple 
mass-concentration relation given by 



„KCDM 



(M,z) 



1 



Mo 
M* 



-0.13 



(18) 



where c V i r = R V i r /r s with r s being the characteris- 
tic scale of an NFW halo, Mo is the halo mass at the 
observed rcdshift z, and M* 



1.5 x 10 13 M m derived 



from simulations [6l|. Thus given a halo mass Mo, we 
can calculate the expected concentration parameter from 
Eq. (fT8"|) . We take a = 1, i.e., z = as the final epoch. 
For Mo = 10 15 M Q , 10 14 M , and 10 13 M , we have 
c vir DM ~ 5.2, 7.0 and 9.5, respectively. For dynami- 
cal dark energy models, the formation and evolution of 
dark matter halos can be different from those in ACDM 
models. Studies show that the differences can largely be 
attributed to the differences of the linear growth factor 
D(z) of the dark matter density perturbations through 
its effect on the halo formation epoch (e.g., [64j][65j]). It 
is found that the concentration parameter c V i r for a dy- 
namical dark energy model can be described well by the 
relation 



ACDM 



D(2 



-oo)/D 



ACDM 



(z -> +oo), (19) 



where D(z — > +oo) and D ACDM (z — > +oo) are the linear 
growth factor evaluated at very high redshift z — > oo for 
the concerned model and for the ACDM model, respec- 
tively ([S^lsHl). Here we adopt this relation to calculate 
the concentration parameter for a particular dark energy 
model with c^f r DM determined by EcLflTS]). The linear 
growth factor is calculated by (e.g., |6q ]) 



D + 2HD - -H 2 tt DM (t)D = 0, 



(20) 



where the initial conditions at very high redshift are 
taken to be D = C x a, and dD/dt = C x da/dt with the 
constant C determined by the normalization condition 
D(z = 0) = 1. With the obtained c V i r , the correspond- 
ing parameter a c / ao is found from Figure 4 of Lu et al. 
H3). We list them in Table H 

Based on Lu et al. [4l| , we use the following methodol- 
ogy to find the initial dark matter mass distribution that 
can give rise to the mass accretion shown in Eq. (|1 T|) . 
For a halo with mass Mo at ao, we assume that M(a) in 
Eq. p7|) is the mass within Mo that is virialized at the 
epoch a. According to the spherical collapse model, the 
linear extrapolated density perturbation averaged over 
the scale corresponding to M(a) should reach the col- 
lapse threshold S c at time a. Then the corresponding 
perturbation at the initial epoch ai can be written as 



Si{M) = S c 



D{ai) 
D[a(M)] 



(21) 



where a(M) is the inverse function of Eq. (fTT)) . repre- 
senting the epoch when the amount of mass M within 



the halo is virialized, and Si(M) is the initial value of 
the density perturbation averaged over the scale within 
which the contained mass is Al. The subscript i denotes 
the quantities at the initial epoch. The initial radius of 
the region with the average density perturbation Si(M) 
is then given by 



n(M) 



3M 



4np D M(ai)[l + 5i(M)] 



1/3 



(22) 



where pdm{ch) is the background dark matter density at 
cii. The initial radial velocity is set by 



Vi(M) 



1 



-H, 



din 5 
din a 



(23) 



where the second term is the peculiar velocity due to the 
mass density perturbation [(37} [HI] . 

We choose the initial rcdshift Zi = 30 (correspondingly 
<H = &o/31), and the initial dark energy perturbations are 
set to be zero. The cosmological parameters are Qdmo = 
0.3, Qdeo = 0.7, and H = 70km/s/Mpc, where Qdmo 
and Qdeo are the current dark matter and dark energy 
density parameters of the universe in unit of the critical 
density. The value of S c is set to be 1.686. 



1.00x10 M,., a c /a =0.8 (F0) A 

c=5.7 

l.O4xlO l4 M a c /a„=0.5 (F6) □ 
c=7.2 

0.82xl0 13 M o 3,730=0.38 (F7) o 

c=9.9 




FIG. 1. Dark matter halo density profiles. The symbols are 
the results from simulations for models of F0 (triangles), F6 
(squares) and F7 (circles), respectively. The lines are the 
corresponding fitted NFW profiles with the concentration pa- 
rameters shown in the plot. 

In Figure [TJ we show the density profiles from our 
1-D simulations for three different halos. The viri- 
alized masses at z = from the simulations are 
M = 1.00 x 10 15 M Q , 1.04 x 10 14 M , and 0.82 x 
10 13 M , respectively, in excellent accordance with Ma = 
10 15 M Q ,10 14 M Q , and 10 13 M , the chosen masses of 
the halos used to set the values of the parameter a c . The 
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TABLE I. Parameters of simulations 





ic 


c\ 


a c /a 




Rvir 


Cvir 










(10 15 M o ) 


(Mpc) 


{_Rvir 1 


FO 


-0.9 


0.5 


0.8 


1.00 


2.6 


5.7 


Fl 


-0.8 


0.5 


0.7 


1.06 


2.6 


5.1 


F2 


-0.7 


0.5 


0.67 


1.08 


2.7 


5.6 


F3 


-0.9 


0.00001 


0.8 


1.10 


2.7 


4.9 


F4 


-0.9 


0.05 


0.8 


1.10 


2.7 


5.2 


F5 


-0.9 


1 


0.8 


1.09 


2.7 


5.3 


F6 


-0.9 


0.5 


0.5 


0.104 


1.2 


7.2 


FT 


-0.9 


0.5 


0.38 


0.0082 


0.52 


9.9 


SQ 


2exp 




0.7 


1.23 


2.8 


5.7 


FQ 


2exp 




0.7 


1.25 


2.8 


5.7 



symbols are the results from the simulations, and the 
lines are the corresponding best-fit NFW profiles with 
the fitted concentration parameters listed in the plot. 
The virial radius R V i r is defined as the radius within 
which the average mass density of the halo is A V i r pDM- 
Here we take A V i r = 337 obtained for the ACDM model 
with n d mo = 0-3 [fUj- Studies have shown that both 5 C 
and A v i r are very insensitive to dark energy models (e.g, 
[iH ) . It can be seen that the simulations with the setting 



described in this section indeed give rise to dark matter 
halos with desired properties. 

To investigate the influence of different physical ingre- 
dients on the behavior of dark energy perturbations, we 
run simulations with different parameters as listed in Ta- 
ble |H F0 to F7 are simulations with the dark energy 
component described by an ideal fluid with the values 
of w and c 2 given in the second and third column. The 
M v i r , R v i r and c V i r are the virial mass, virial radius and 
the fitted NFW concentration parameter measured from 
simulations. The last two models are for the quintessence 
dark energy with the double exponential potential given 
by Eq. (fl"4)) . The 'SQ' run calculates the dark energy per- 
turbations by considering the scalar field evolution shown 
in Eq. (JSJ), and the 'FQ' run uses the fluid approach with 
the corresponding w and c 2 to that of the quintessence 
model (see Eq. © and Eq. (fl"3|) ). 



IV. RESULTS 

In this section, the results on dark energy perturba- 
tions and their dependence on different physical param- 
eters are presented. The model F0 is taken to be the 
reference model. 



A. Reference model 

We first show the temporal and spatial behaviors of 
dark energy perturbations for F0. 



The left panel of Figure [2] shows the redshift evolution 
of the dark matter density contrast 5dm and the dark 
energy perturbation Sde for the most inner radial-mesh 
bin r = 0.1 Mpc. It can be seen that for the central 
part, the dark matter virialization happens at z ~ 4. 
The value of 5dm reaches ~ 10 5 at z = 0. For Sde, 
it is always positive and increases smoothly without a 
sharp change corresponding to the virialization of dark 
matter at z ~ 4. Note that for F0, c 2 = 0.5, which gives 
rise to a resistance to the perturbation growth for the 
dark energy component and explains its different time 
evolution behavior in comparison with that of the dark 
matter. At z = 0, we have Sde ~ 10~ 5 , which is about 
ten orders of magnitude smaller than 5 dm- 

The right panel of Figure [2] shows the spatial profiles 
for 5dm and Sde, respectively, at different z. At z > 5, 
the halo is in its early formation stage, and no dark mat- 
ter virializations occur. As the evolution proceeds, the 
virializations start from the central region and continu- 
ally extend to larger regions. At z = 1.2, the virialized 
part reaches r ~ 0.8 Mpc as indicated by the change of 
the 5dm profile. At z = 0, the virial region extends to 
r « 2.6 Mpc, fully consistent with the virial radius de- 
fined by the average density contrast A„j r ss 337 with 
respect to the average matter density of the universe for 
Obmo = 0.3. For the dark energy perturbation Sde, its 
amplitude grows with the decrease of the redshift, and 
the profile gets steeper. As discussed above, unlike the 
behavior of 5dm from which the virialized and unvirial- 
ized regions can be easily identified, 5de is rather smooth 
and there is not an apparent feature at the dark matter 
virialization boundary. Furthermore, the profile of 5de 
is much shallower than that of the dark matter. From 
10 Mpc, 5de decreases from 10 -5 to 



r = 0.2 Mpc to r = 
10~ 6 , whereas for 



'dm 



, it changes from ~ 10 to 0(1). 



B. Dependences of Sde on w and cj 

In Figure [31 we show the dependence of 5de on the 
equation of state parameter w of the dark energy com- 
ponent. The time evolution of 5de for the inner- most 
radial-mesh bin is presented for w = —0.9 (F0), w = —0.8 
(Fl) and w = —0.7 (F2), respectively. The parameter c 2 
is fixed to be c 2 = 0.5. The ^-dependence is clearly seen. 
For w = —0.7, the dark energy perturbation 5de is about 
three times as large as that with w = —0.9. 

Figure [J] shows the effects of c 2 on dark energy per- 
turbations with c 2 = 1.0 (F5), 0.5 (F0), 0.05 (F4) 
and 0.00001 (F3), respectively. A sensitive dependence 
of Sde on c 2 is apparent. For the three cases with 
c 2 > 0.05, the dark energy perturbation amplitude in- 
creases with time smoothly without characteristic fea- 
tures corresponding to the virialization of the dark mat- 
ter component. This is because the large sound speed of 
the dark energy can coordinate its behavior over a large 
range quickly in comparison with that of the dark matter. 
As expected, the amplitude of Sde is larger for smaller c 2 . 
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FIG. 2. Left panel: The redshift evolution of the dark energy perturbation Sde (solid line)for the radial-mesh bin r — 0.1 Mpc 
is shown for reference model FO (M = 1O 15 M0, w = —0.9, (? a = 0.5). The corresponding dark matter density perturbation 5dm 

(dotted lines) are shown for model 



is also shown (dotted line). 
FO at z = 0, z = 1.2 and z 



Right panel: The spatial profiles of Sde (solid lines) and 5dm 
- 5.4, respectively. 




FIG. 3. The redshift evolution of Sde for the radial-mesh bin 
r = 0.1 Mpc is shown for F0 (solid line), Fl (dotted line) and 
F2 (dot-dashed line), respectively. 

For c 2 = 0.05, Sde is about an order of magnitude higher 
than that of c 2 = 0.5. For a very small c 2 = 0.00001, we 
can see that at z ~ 4, there is a weak feature indicat- 
ing the virialization of the dark matter component. The 
amplitude of Sde in this case reaches 0(1) at z = in 
the central region. In outer parts of the halo and at high 
redshifts, we still have Sde < 1- Thus our linear analy- 
ses for dark energy perturbations are still approximately 



FIG. 4. The c^-dependence of Sde- The redshift evolution of 
Sde for the inner-most bin at r = 0.1 Mpc is shown for F3 
(dotted line), F4 (dashed line), F5 (dot-dashed line) and F0 
(solid line). 

valid for c 2 = 0.00001. For even lower c 2 , we expect that 
the dark energy perturbation would follow that of the 
dark matter more closely and can reach an amplitude of 
Sde > 1- In such cases, our formulation presented in 
Sec. II for linear dark energy perturbations are not ap- 
plicable, and nonlinear dark energy perturbations must 
be considered carefully. Moreover, the feedback effects of 
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dark energy perturbations on the formation of the dark 
matter halo may also need to be taken into account (e.g., 

Hi). 



Dependence of Sde on properties of dark 
matter halos 



Here we study the dark energy perturbations induced 
by different dark matter halos. In the left panel of Fig- 
ure [5j the time evolution of Sde for the inner-most bin 
is shown for M = 1.00 x 10 15 M (F0), 1.04 x 10 14 M Q 
(F6), and 0.82 x 10 13 M© (F7), respectively. It is seen 
that more massive halos induce larger dark energy per- 
turbations. At z = 0, Sde of F0 is about 20 times larger 
than that of model F7. 

The right panel of Figure [5] shows the spatial profiles 
of Sde for the three halos at z = 5.4 and 0, respectively. 
It is noted that for smaller halos, the dark matter den- 
sity profile is more concentrated, and occupies smaller 
regions. Consequently, the profile for Sde is steeper and 
extends less for less massive halos. On the other hand, in 
all the three cases, the dark energy perturbations extend 
to regions much larger than the virial radii of the halos. 
The somewhat different line shapes in the left panel of 
Figure [5] are related to the different virialization epochs 
for the central region for different models. The less mas- 
sive halo virializes earlier. 

We further analyze the approximate relation between 
Sde and the halo properties. On halo scales, we have 
k/a 3> H. From the dark energy perturbation equations 
Eq. (0-([6]), we thus have approximately 



f 



-S DE + $ ~ 0, S DE 



1 



w 



(24) 



In Figure [BJ we show Sde/S* at z = for a number 
of models with different mass of dark matter halos, dif- 
ferent w and c 2 for the dark energy component, where 
5* = — (1 + w)$/c 2 . It is seen that for all the mod- 
els, \5de/S* — 1| < 1% within the virial radius and 
\Sde/S* — 1| < 4% for r < 5R V i r . This demonstrates that 
the dark energy perturbation traces the gravitational po- 
tential of dark matter halos very well over a broad range 
of model parameters, which may provide us an efficient 
way to estimate the power spectrum of dark energy per- 
turbations in the regime of nonlinear structure forma- 
tion. A similar tight correlation between Sde and <3? is 
also shown in [25[ for the extended quintessence model. 



D. Correspondence between the fluid approach 
and the scalar field model 



As discussed in Sec. II, the correspondence between 
a scalar field dark energy model and its fluid description 
can be found by specifying suitable w and c 2 that depend 
on the dynamical evolution of the scalar field. Such cor- 
respondences have been applied extensively in analyzing 



the effect of the background dark energy on the structure 
formation and the behavior of the dark energy perturba- 
tion in the linear regime of the structure formation (e.g., 
14|[z3|)- Here we test their validity in studying dark 
energy perturbations induced by the nonlinear structure 
formation. We consider the scalar field model with a dou- 
ble exponential potential given in Eq. (| 14[) with a = 6 
and ft = 0.1 , in accord with our previous studies 31 



In the simulation runs SQ and FQ, we calculate the dark 
energy perturbation from the scalar field perturbation 
given in Eq. Q and from the fluid approach with w and 
c 2 s changing with time according the dynamical evolution 
of the scalar field, respectively. 

The time evolution of Sde for the inner- most bin is 
shown in the upper panel of Figure [7] for SQ and FQ, re- 
spectively. The ratio of the two r$ = Sde{SQ)/Sde(FQ) 
is shown in the lower panel. The corresponding spatial 
behaviors at z = 5.3,1.2 and z = are shown in Fig- 
ure [SJ It can be seen that the results from SQ and 
FQ match very well. The relative large scatter with 
rs ~ 20% at low redshift for the inner-most bin are due 
to the isotropized virialization that induces random dis- 
persions for the gravitational potential somewhat differ- 
ently in different model runs. These comparisons demon- 
strate that the fluid approach can be applied to study the 
dark energy perturbation behavior in the whole regime 
of structure formation, from linear to nonlinear stages. 



V. SUMMARY AND DISCUSSION 

With two-component 1-D numerical simulations tak- 
ing into account the fast and slow growth of dark matter 
halos, wc analyze the behavior of dark energy perturba- 
tions induced by the formation of spherical dark matter 
halos. By comparing the results calculated directly from 
the scalar field evolution and that from the fluid descrip- 
tion for a quintessence dark energy model with a double- 
exponential potential, we show that the fluid approach 
can be used to analyze the dark energy behavior very well 
even in the nonlinear stage of structure formation. In the 
fluid treatment for dark energy, the equation of state w 
and the sound speed c s are the two important parameters 
that affect the dark energy perturbations significantly. 
Our studies find that in general, the dark energy pertur- 
bation arising from the formation of a dark matter halo 
has a much more extended profile than that of the dark 
matter halo except for the case c s where the dark 
energy component follows the dark matter closely. A re- 
lation of Sde ~ [—(1 + w)/c 2 ]§ with an accuracy about 
1 — 2% within the virial radius of a halo is revealed, which 
provides us a potential means to estimate the power spec- 
trum of dark energy perturbations from that of the dark 
matter potential field. From w = —0.9 to w = —0.7, 
Sde increases about three times from ~ 2.5 x 10~ 5 to 
7.5 x 10~ 5 for c 2 = 0.5. On the other hand, varying c 2 
from 0.5 to 0.05, the dark energy perturbation increases 



by an order of magnitude from ~ 10 



to 



io- 



Our 
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FIG. 5. The left panel shows the dependence of Sde on the mass of the dark matter halo. The evolution of Sde at r = 0.1 Mpc 
is shown for F0, F6 and F7, respectively. The right panel shows the spatial distribution of Sde for the three models at z — 5.4 
(the lower set) and z = (the upper set), respectively. 



c s 2 =0.05 (F4) 
c, 2 =0.5 (F0) 
c s 2=1.0 (F5) 
c s 2 =0.5 (F2) 
c 5 2 =0.5 (F7) 



FIG. 6. The validity of the relation given by Eq. (|24|) between 
See and $ is shown for different models. 6* = —(1 + w)&/Cg. 



analyses show that our simulation treatment and the lin- 
ear evolution of dark energy perturbations are valid for 
Cg down to Cg ~ 0.00001. For even smaller c% more accu- 
rate analyses taking into account nonlinear dark energy 
perturbations and their feedback effects on the formation 
of dark matter halos are needed (e.g.,[3!|). The depen- 
dence of the dark energy perturbation on the mass of the 
halo is also analyzed. It is shown that Sde increases by 
- 20 times from M « 10 13 M Q to M « 10 15 M . 

The numerical analyses presented in this paper are 
done for dark energy models with w > — 1. Our sim- 
ulations show that for such models, the dark energy per- 
turbation induced by the formation of dark matter halos 
has a clustering behavior with Sde > during the entire 
evolutionary process. This can also be seen clearly from 
the approximate relation Sde ~ [— (1 + w )/ci]® (note 




:o 25 ju 



FIG. 7. The redshift dependence of Sde for the inner-most 
bin is shown in the upper panel for SQ and FQ, respectively. 
In the lower panel, the ratio rs of the two is shown. 



$ < in dark matter halo regions). On the other hand, 
for w < —1, dark energy voids are expected from this 
relation although wc do not study these cases explicitly 
in our simulations. These results are in good agreement 
with studies shown in, e.g., 34 1 35|. However, there are 



other analyses that point to the existence of dark en- 
ergy void with Sde < during the quasi-linear stages of 
dark matter halo formation even for models with w > — 1 
(e.g., [29] 30] [31.]). It is noticed that all the latter studies 
are performed using the Lemaitre-Tolman-Bondi (LTB) 
metric that is in the synchronous gauge. On the other 
hand, our simulation analyses presented here adopt the 
Newtonian gauge. It is well known that cosmological 
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FIG. 8. The spatial behaviors of Sde for SQ and FQ are 
shown in the upper panel for z — 5.3, 1.2 and z = 0, respec- 
tively. The corresponding are shown in the lower panel. 



72114c 



b~DE(Syn) = 5 DE (New) - a 



Pde 
Pde 



(25) 



where a(t, x) is related to the coordinate transforma- 
tion for the time component associated with the two 
gauges, t(Syn) = t(New) + a, and da/dt — $. We have 
Pde I Pde = — 3(1 + w)H. Further with the approximate 
relation between 5de(Ncw) and $, we then obtain 



WE 



(New) 



dt. 



ode (Syn) « 6 be (New) - 3(1 + w)Hc 2 s 

(26) 

Thus if the second term is larger than the first term, 
a dark energy void with dDE(Syn) < can occur in 
the LTB synchronous gauge even for w > — 1 with 
SoE(New) > 0. In other words, our results about the 
dark energy clustering for w > —1 presented in this pa- 
per do not conflict with those showing the existence of 
dark energy voids for the same models. Instead, the dif- 
ferences merely reflect the differences of the specific gauge 
used in different analyses. 



energy density perturbations are gauge dependent (e.g., 



71 



[43|[72j). Concerning the dark energy density pertur- 
bation in the LTB synchronous gauge SuE(Syn) and in 
the Newtonian gauge SoE(New), we have the relation 
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